%Modified 2/19/11 to make a track-by-track graph of cell length vs. time,
%for Thomas Murooka

function [varargout] = GraphAnalysis(min_inst_vel, cyc_time, varargin)

% [varargout] = GraphAnalysis(infile)
% Inputs:
%  infile - Can either be a single string of a filename to plot, or a cell
%       array of filenames.

warning off all;
num_files = length(varargin);
filelist = varargin;

% Create the figures for plotting
distrib_fig = figure; hold on;
axis square;
set(gcf, 'color', 'white');
set(gca, 'color', 'white');
set(gcf, 'InvertHardCopy', 'off');
distrib_fig_two = figure; hold on;
axis square;
set(gcf, 'color', 'white');
set(gca, 'color', 'white');
set(gcf, 'InvertHardCopy', 'off');
cell_mot_chg_fig = figure; hold on;
axis square;
set(gcf, 'color', 'white');
set(gca, 'color', 'white');
set(gcf, 'InvertHardCopy', 'off');
scatter_fig = figure; hold on;
axis square;
set(gcf, 'color', 'white');
set(gca, 'color', 'white');
set(gcf, 'InvertHardCopy', 'off');
cell_fft_fig = figure; hold on;
axis square;
set(gcf, 'color', 'white');
set(gca, 'color', 'white');
set(gcf, 'InvertHardCopy', 'off');
% scatter_fig_two = figure; hold on;
% set(gcf, 'color', 'white');
% set(gca, 'color', 'white');
% set(gcf, 'InvertHardCopy', 'off');
contour_plot = figure; hold on;
axis square;
set(gcf, 'color', 'white');
set(gca, 'color', 'white');
set(gcf, 'InvertHardCopy', 'off');
change_hist_fig = figure; hold on;
axis square;
set(gcf, 'color', 'white');
set(gca, 'color', 'white');
set(gcf, 'InvertHardCopy', 'off');
sample_fig = figure; hold on;
axis square;
set(gcf, 'color', 'white');
set(gca, 'color', 'white');
set(gcf, 'InvertHardCopy', 'off');
fft_fig = figure; hold on;
axis square;
set(gcf, 'color', 'white');
set(gca, 'color', 'white');
set(gcf, 'InvertHardCopy', 'off');
anglechg_fig = figure; 
set(gcf, 'color', 'white');
set(gca, 'color', 'white');
set(gcf, 'InvertHardCopy', 'off');
% plot added by Francesco 2/19/2011
length_tracks = figure; hold on;
set(gcf, 'color', 'white');
set(gca, 'color', 'white');
% plot added by Francesco 6/29/2011
vel_tracks = figure; hold on;
set(gcf, 'color', 'white');
set(gca, 'color', 'white');
% plot added by Francesco 6/29/2011
length_vel = figure; hold on;
set(gcf, 'color', 'white');
set(gca, 'color', 'white');

count = 0;
colors = {'r', 'g', 'b', 'm', 'k', 'y', 'c'};
cont_colors = {'r', 'g', 'b', 'c', 'm', 'k', 'y'};

hue0 = [0, 0, 1];  % palette added by Francesco Marangoni 6/29/2011 - blue 
hue1 = [0, 0.3, 1];
hue2 = [0, 0.6, 1];
hue3 = [0, 1, 1];  % teal
hue4 = [0, 1, 0.6]; 
hue5 = [0, 1, 0.3];
hue6 = [0, 1, 0];   % green
hue7 = [0.3, 1, 0];
hue8 = [0.6, 1, 0];
hue9 = [1, 1, 0];   % yellow
hue10 = [1, 0.6, 0];
hue11 = [1, 0.3, 0];
hue12 = [1, 0, 0];   % red

wh = waitbar(0, 'Graphing file(s)...');

for cur_file = 1:num_files,
    clear infile;

    infile = filelist{cur_file};
    
    fid = fopen(infile);
    
    % init variables
    all_data = {};
    num_tracks = 1;
    numlines = 1;
    data = [];
    
    % ignore first line
    tline = fgetl(fid);
    
    % read data out of file
    while ~feof(fid),        
        % get the whole line
        tline = fgetl(fid);
        
        if (strcmp(upper(tline(1:3)), '---')),
            % so, we are at the end of a cell's track
            all_data{num_tracks} = data;
            numlines = 1;
            data = [];
            lineofdata = {};
            num_tracks = num_tracks + 1;
            
        else
            % read in lines of data for a given track
            datacount = 1;
            lineofdata = {};
            % read in the chunks of the line
            while (length(tline) > 0),
                [token, tline] = strtok(tline, char(9));
                lineofdata{datacount} = token;
                datacount = datacount+1;
            end
            
            % transfer the chunks into numbers
            for i = 1:24,  %modified 2/19/11 to handle column 23 (delta time) and 24 (cell length)
                try
                    data(numlines, i) = str2num(lineofdata{i});
                catch
                    data(numlines, i) = nan;
                end
            end
            numlines = numlines + 1;
        end
    end
    
    % Scheme of length vs. time by tracks in 2D Francesco Marangoni 2/19/2011
    
    
    
    grey = [0.9 0.9 0.9];
    max_length = input('Insert the max length to plot on Y (microns)');
    max_vel = input ('Insert the max velocity to plot on Y (microns/min)');
    for i = 1:length(all_data),
        time = [];
        cell_length = [];
        mod_cell_length = [];
        cell_vel = [];
        mod_cell_vel = [];
        signhue = zeros(1,3);
        temp = all_data{i};
        time = [time; temp(find(~isnan(temp(:,22))),22)];
        cell_length = [cell_length; temp(find(~isnan(temp(:,24))),24)];
        
        mod_cell_length = [mod_cell_length; temp(find(~isnan(temp(:,24))),24)];
        
        mod_cell_length(:) = (2*i+(mod_cell_length(:)./max_length)-1);
        cell_vel = [cell_vel; temp(find(~isnan(temp(:,10))),10)];
        mod_cell_vel = [mod_cell_vel; temp(find(~isnan(temp(:,10))),10)];
        mod_cell_vel(:) = (2*i+(mod_cell_vel(:)./max_vel)-1);
        cell_vel = [NaN; cell_vel];  %adds NaN at the beginning of the vector
        mod_cell_vel = [NaN; mod_cell_vel];  %adds NaN at the beginning of the vector
        el_time = zeros(length(time)); 
        ranking = zeros(length(time));
        ranking(:) = 2*i; 
        for j = 1:length(time),
            el_time(j) = (time(j)-1)*cyc_time/60;  %plots actual track start
        end
        
        figure(length_tracks);% opens length vs. time figure
        rectangle('Position', [el_time(1), 2*i-1, (el_time(length(el_time))-el_time(1)), 1], 'EdgeColor', grey, 'FaceColor', grey);
        for step = 2:(length(el_time)),  %jumps first velocity value because it is NaN
            if cell_vel(step)<=2,  % PALETTE ASSIGNMENT for length vs. time
                signhue = hue0;
            elseif (cell_vel(step)>2 & cell_vel(step)<=4),
                signhue = hue1;
            elseif (cell_vel(step)>4 & cell_vel(step)<=6),
                signhue = hue2;
            elseif (cell_vel(step)>6 & cell_vel(step)<=8),
                signhue = hue3;
            elseif (cell_vel(step)>8 & cell_vel(step)<=10),
                signhue = hue4;
            elseif (cell_vel(step)>10 & cell_vel(step)<=12),
                signhue = hue5;
            elseif (cell_vel(step)>12 & cell_vel(step)<=14),
                signhue = hue6;
            elseif (cell_vel(step)>14 & cell_vel(step)<=16),
                signhue = hue7;
            elseif (cell_vel(step)>16 & cell_vel(step)<=18),
                signhue = hue8;
            elseif (cell_vel(step)>18 & cell_vel(step)<=20),
                signhue = hue9;
            elseif (cell_vel(step)>20 & cell_vel(step)<=22),
                signhue = hue10;
            elseif (cell_vel(step)>22 & cell_vel(step)<=24),
                signhue = hue11;
            else
                signhue = hue12;
            end
            set(gca, 'YTickLabel', [0; max_length]);
           
            plot(el_time(step), mod_cell_length(step), '.', 'Color', signhue, 'Linewidth', 5);
        end
        
        for step = 1:(length(el_time)-1),
            plot(el_time(step:(step+1)), mod_cell_length(step:(step+1)), 'Color', 'k', 'Linewidth', 1);
        end
        
        figure(vel_tracks);%opens velocity vs time figure
        rectangle('Position', [el_time(1), 2*i-1, (el_time(length(el_time))-el_time(1)), 1], 'EdgeColor', 'w', 'FaceColor', 'w');

        for step = 2:(length(el_time)-1),%jumps first value because velocity is NaN
            plot(el_time(step:(step+1)), mod_cell_vel(step:(step+1)), 'Color', 'k', 'Linewidth', 1);
        end

        for step = 2:(length(el_time)),  %jumps first value because velocity is NaN
            if cell_length(step)<=20,  % PALETTE ASSIGNMENT for velocity vs. time
                signhue = hue0;
            elseif (cell_length(step)>20 & cell_length(step)<=21),
                signhue = hue1;
            elseif (cell_length(step)>21 & cell_length(step)<=22),
                signhue = hue2;
            elseif (cell_length(step)>22 & cell_length(step)<=23),
                signhue = hue3;
            elseif (cell_length(step)>23 & cell_length(step)<=24),
                signhue = hue4;
            elseif (cell_length(step)>24 & cell_length(step)<=25),
                signhue = hue5;
            elseif (cell_length(step)>25 & cell_length(step)<=26),
                signhue = hue6;
            elseif (cell_length(step)>26 & cell_length(step)<=27),
                signhue = hue7;
            elseif (cell_length(step)>27 & cell_length(step)<=28),
                signhue = hue8;
            elseif (cell_length(step)>28 & cell_length(step)<=29),
                signhue = hue9;
            elseif (cell_length(step)>29 & cell_length(step)<=30),
                signhue = hue10;
            elseif (cell_length(step)>31 & cell_length(step)<=300),
                signhue = hue11;
            else
                signhue = hue12;
            end
            set(gca, 'YTickLabel', [0; max_length]);
            plot(el_time(step), mod_cell_vel(step), '.', 'Color', signhue, 'Linewidth', 5);
        end
        
      
        
        
        
    end
    figure(length_tracks);% opens length vs. time figure
    set(gca,'YTick',[1:1:(2*length(all_data))]);
    %axis ([0, 60, 0, 2*(length(all_data)+1)]);
    xlabel('Time (min)');
    ylabel('Cell length (\mum)');
    title ('Schematic view Cell Length vs. Time in 2D');
    
    figure(vel_tracks);%opens velocity vs time figure
    set(gca,'YTick',[1:1:(2*length(all_data))]);
    %axis ([0, 60, 0, 2*(length(all_data)+1)]);
    xlabel('Time (min)');
    ylabel('Instantaneous velocity (\mum/min)');
    title ('Schematic view Cell Velocity vs. Time in 2D');
  
   % Scatterplot Inst Vel 3D vs. Cell length  Francesco Marangoni 6/29/11
    figure(length_vel);
    x1 = [];
    y1 = [];
    for i = 1:length(all_data),  % puts together velocity with length
        temp = all_data{i};
        y1 = [y1; temp(find(~isnan(temp(:,24))),24)];
        y1(1) = [];  %removes first element of vector
        x1 = [x1; temp(find(~isnan(temp(:,10))),10)];
    end
    
    plot (x1,y1,'bo','MarkerFaceColor','b','MarkerSize',4);
    xlabel('Inst. 3D velocity (\mum/min)');
    ylabel('Cell length (\mum)');
    title ('Instantaneous 3D velocity vs. Cell length');
    grid on;
    
    
    % 3D velocity histogram
    x = 0:2:40;
    y = [];
    for i = 1:length(all_data),
        temp = all_data{i};
        y = [y; temp(find(~isnan(temp(:,10))),10)];
    end
%     [phat, pci] = gamfit(y);
%     [m,v] = gamstat(phat(1), phat(2));
%     [lambdahat,lambdaci] = poissfit(y);
%     [weibp, weibci] = weibfit(y);
%     fprintf('Distribution %d: %f (%f - %f); %f (%f - %f)\n', cur_file, phat(1), pci(1,:), phat(2), pci(2,:));
%     fprintf('\tMean - %f; Variance - %f\n', m, v);
%     fprintf('\tLambda - %f (%f - %f)\n', lambdahat, lambdaci);
%     fprintf('\tWeibull - %f (%f - %f); %f (%f - %f)\n', weibp(1), weibci(1,:), weibp(2), weibci(2,:));
    for_ttest{cur_file} = y;
    n = hist(y,x);
    distrib_all(cur_file, :) = n/length(y);
        
    % Scatterplot instantaneous velocity x versus instantaneous velocity y
    figure(scatter_fig);
    x = [];
    y = [];
    for i = 1:length(all_data),
        temp = all_data{i};
        y = [y; temp(find(~isnan(temp(:,12))),12)];
        x = [x; temp(find(~isnan(temp(:,11))),11)];
    end
    plot(x,y, '.', 'Color', colors{mod(cur_file-1,7)+1});
    xlabel('Inst. velocity in x (\mum/s)');
    ylabel('Inst. velocity in y (\mum/s)');
    %title(sprintf('Instantaneous velocity in x vs. Instantaneous velocity in y - %s', strrep(strtok(infile,'.'), '_', '\_')));
    title('Instantaneous velocity in x vs. Instantaneous velocity in y');
    axis([-40 40 -40 40]);
    grid on;
    
%     figure(scatter_fig_two);
%    [xmu, xsig, xmuci, xsigci] = normfit(x);
%   [ymu, ysig, ymuci, ysigci] = normfit(y);
    xplot = -40:2:40;
    yplot = -40:2:40;
    clear zz;
    clear z;
%     for g = 1:length(xplot),
%         for h = 1:length(yplot),
%             z(g,h) = normpdf(xplot(g), xmu, xsig) * normpdf(yplot(h), ymu, ysig);
%         end
%     end
    for g = 1:length(xplot)-1,
        for h = 1:length(yplot)-1,
            tempx = logical(x >= xplot(g)) & logical(x <= xplot(g+1));
            tempy = logical(y >= yplot(h)) & logical(y <= yplot(h+1));
            zz(g,h) = sum(ismember(find(tempx), find(tempy)));
        end
    end
%     [tx, ty] = meshgrid(xplot, yplot);
%     contour3(tx, ty, z);
%     xlabel('Inst. velocity in x (\mum/s)');
%     ylabel('Inst. velocity in y (\mum/s)');
%     %title(sprintf('Instantaneous velocity in x vs. Instantaneous velocity in y - %s', strrep(strtok(infile,'.'), '_', '\_')));
%     title('Instantaneous velocity in x vs. Instantaneous velocity in y -- Gaussian Fit');
%     axis([-40 40 -40 40]);
%     grid on;
    
    %contour plot
    figure(contour_plot);
    [tx, ty] = meshgrid(-39:2:39, -39:2:39);
    contour3(tx, ty, zz, 20, cont_colors{mod(cur_file-1, 7)+1});
    xlabel('Inst. velocity in x (\mum/s)');
    ylabel('Inst. velocity in y (\mum/s)');
    %title(sprintf('Instantaneous velocity in x vs. Instantaneous velocity in y - %s', strrep(strtok(infile,'.'), '_', '\_')));
    title('Instantaneous velocity in x vs. Instantaneous velocity in y');
    axis([-40 40 -40 40]);
    grid on;
    
    
    %2D vel hist data generation
    hx = 0:2:40;
    hist_data_2d{cur_file} = sqrt(x.^2 + y.^2);
    n = hist(sqrt(x.^2 + y.^2), hx);
    distrib_2d(cur_file, :) = n/length(y);    
    
    %Absolute Angle Change histogram (in progress)
    %figure;
    angle = [];
    for i = 1:length(all_data),
        temp = all_data{i};
        for j = 1:size(temp,1),
            if (((temp(j,7)) == 0) || (temp(j,8)) == 0),
                angle = [angle; 0];
            else
                temp_t = atan(temp(j,8)/temp(j,7));
                temp_t = (temp_t/pi)*180;
                if (~isnan(temp_t)),
                    angle = [angle; temp_t];
                end
            end
        end
        
    end
    x =-180:10:180;
    n = hist(angle, x);
    %bar(n)
    %title('Absolute Angle Histogram');
    %output absolute angle array
    %filename = strcat(infile,'absangle.txt');
    %dlmwrite(filename, n, '\t');
    
    
    % Angle Change histogram
    figure(change_hist_fig);
    x = -180:10:180;
    y = [];
    for i = 1:length(all_data),
        temp = all_data{i};
        for j = 1:size(temp,1),
            if ~isnan(temp(j,20)),
                y = [y; temp(j,20)];
            end
        end
    end
    n = hist(y, x);
    all_change(cur_file, :) = n/length(y);
    
 


    
    %polar plot of angle change
    figure(anglechg_fig);
    dir = (x * pi/180);
    polar(dir, n, colors{mod(cur_file-1, 7)+1})
    hold on;
    title('Angle Change');
    
    % Representative Tracks (need to plot w/all starting from 0)
    figure(sample_fig);
    hold on;
    x = [];
    y = [];
    axis ([-100 100 -100 100]);
    set (gca, 'FontSize', 18, 'LineWidth', 2);
    for i = 1:length(all_data),
        temp = all_data{i};
        y = [0; temp(:,7)];
        x = [0; temp(:,8)];
        track = zeros(size(temp,1), 2);
        for j = 2:size(temp,1),
            track(j,1) = track(j-1,1)+temp(j,8);
            track(j,2) = track(j-1,2)+temp(j,7);
        end    
        if (num_files == 1),
            plot(track(:,1), track(:,2), 'Color', colors{mod(count, 7)+1}, 'LineWidth', 1);
        else
            plot(track(:,1), track(:,2), 'Color', colors{mod(cur_file-1, 7)+1}, 'LineWidth', 1);
        end
        count = count + 1;

        xlabel('\Deltax (\mum)');
        ylabel('\Deltay (\mum)');
        title('Cell Tracks');


    end
    [legend_h,object_h,plot_h,text_strings] = legend(strrep(filelist,'_','\_'), 'Location', 'SouthOutside');
    for i = 1:length(filelist),
        set(object_h(i), 'FontSize', 12);
    end
    for i = 1:length(plot_h),
        set(plot_h(i), 'Color', colors{mod(i-1, 7)+1});
    end
    
    %Calculate fft for each cell
    figure(cell_fft_fig);
    max = 0;
    for i = 1:length(all_data),
        if (max<size(all_data{i},1)),
            max = size(all_data{i},1);
        end
    end
    max = ceil(max/2)*2;
    count = 0;
    clear fft_cell_array fft_cell_toplot_array fft_array fft_toplot_array;
    for i = 1:length(all_data),
        if (size(all_data{i},1) >= 10),
            count = count + 1;
            temp = all_data{i};
            vel_temp = sum(temp(2:end,11:12).^2,2)';
            fft_array(count, :) = zeros(1, max);
            fft_array(count, 1:length(vel_temp)) = detrend(vel_temp);
        end
    end
    for i = 1:size(fft_array,1),
        Y = fft(fft_array(i,:), max);
        Pyy = abs(Y)/max;
        fft_toplot_array(i,:) = Pyy;
    end    
    % THIS ASSUMES YOUR CYCLE TIME IS 15s (so, a minute requires *4)
    x = (0:max/2)/max*4;
    plot(x,fft_toplot_array(:,1:(max/2)+1), 'Color', colors{mod(cur_file-1, 7)+1}, 'LineWidth', 1.25);   
    xlabel('Frequency (cycles/min)');
    ylabel('Power');
    title('Power Spectrum of Each Cell''s Velocity Change over Time');
    
    %output single cell fft
    k = findstr(infile,'/');
    if isempty(k),
        infile = infile(1:end-12);
    else
        infile = infile(k(end)+1:end-12);
    end
    filename = strcat(infile,'sglcell_fft.txt');
    dlmwrite(filename, x, 'delimiter','\t');
    dlmwrite(filename, fft_toplot_array(:,1:(max/2)+1), '-append', 'delimiter','\t')

    
    
    % Calculate fft of inst vel x and y for whole files 
    figure(fft_fig);
    max = 0;
    for i = 1:length(all_data),
        if (max<size(all_data{i},1)),
            max = size(all_data{i},1);
        end
    end
    max = ceil(max/2)*2;
    count = 0;
    clear fft_array fft_toplot_array mean_fft_toplot_array;
    for i = 1:length(all_data),
        if (size(all_data{i},1) >= 10),
            count = count + 1;
            temp = all_data{i};
            vel_temp = sum(temp(2:end,11:12).^2,2)';
            fft_array(count, :) = zeros(1, max);
            fft_array(count, 1:length(vel_temp)) = detrend(vel_temp);
        end
    end
    for i = 1:size(fft_array,1),
        Y = fft(fft_array(i,:), max);
        Pyy = abs(Y)/max;
        fft_toplot_array(i,:) = Pyy;
    end
    mean_fft_toplot_array = mean(fft_toplot_array,1);
    
    
    % THIS ASSUMES YOUR CYCLE TIME IS 15s (so, a minute requires *4)
    x = (0:max/2)/max*4;
    plot(x,mean_fft_toplot_array(1:(max/2)+1), 'Color', colors{mod(cur_file-1, 7)+1}, 'LineWidth', 1.25);   
    xlabel('Frequency (cycles/min)');
    ylabel('Power');
    title('Average Power Spectrum of Velocity Change over Time');
    
    fclose(fid);
    
    [outfile, rem] = strtok(filelist{cur_file}, '.');
    outfile = strcat(outfile, '_fft.txt');
    fid = fopen(outfile, 'w');
    fprintf(fid, 'Freq (cyc/min)\tPower\n');
    for i = 1:length(x)
        fprintf(fid, '%5.2f\t%5.2f\n',x(i), mean_fft_toplot_array(i));
    end
    fclose(fid);
    
    waitbar(cur_file/num_files);
end
close(wh);

%distribution of 3DIV 
figure(distrib_fig);
x = 0:2:40;
set (gca, 'FontSize', 18);
bar(x, distrib_all', 1);
plot(x, distrib_all');
xlabel('Inst. 3D velocity (\mum/s)')
ylabel('Relative Frequency')
%title(sprintf('Distribution of 3D Instantaneous Velocities - %s', strrep(strtok(infile,'.'), '_', '\_')));
title('Distribution of 3D Instantaneous Velocities');
axis([-0.5 40.5 0 0.35])
%if num_files == 2,
%    p = ranksum(for_ttest{1}, for_ttest{2}, 0.05);
%    fprintf('Rank Sum Test for 3D Velocities -- p <= %d\n', p);
%end
legend(strrep(filelist,'_','\_'),'Location','SouthOutside');

figure(distrib_fig_two);
x = 0:2:40;
bar(x, distrib_2d', 1);
xlabel('Inst. 2D velocity (\mum/s)')
ylabel('Relative Frequency')
%title(sprintf('Distribution of 3D Instantaneous Velocities - %s', strrep(strtok(infile,'.'), '_', '\_')));
title('Distribution of 2D Instantaneous Velocities');
axis([-0.5 40.5 0 0.35])
%if num_files == 2,
%    p = ranksum(hist_data_2d{1}, hist_data_2d{2}, 0.05);
%    fprintf('Rank Sum Test for 2D Velocities -- p <= %d\n', p);
%end
legend(strrep(filelist,'_','\_'), 'Location','SouthOutside');

figure(scatter_fig);
legend(strrep(filelist,'_','\_'));

figure(fft_fig);
legend(strrep(filelist,'_','\_'));

figure(change_hist_fig);
x =-180:10:180;
bar(x, all_change', 1);
xlabel('Change in Cell Direction (\circ)');
ylabel('Relative Frequency');
%title(sprintf('Change in Cell Direction - %s', strrep(strtok(infile,'.'), '_', '\_')));
title('Change in Cell Direction');
axis([-180.5 180.5 0 0.1]);
legend(strrep(filelist,'_','\_'));

%figure(sample_fig);
%legend(filelist);


% %Migrating Cells

all_times = [];
for i = 1:length(all_data),
    temp = all_data{i};
    all_times = [all_times; temp(2:end,22)];    
end
times = unique(all_times(find(~isnan(all_times))));
pop_vel = cell(size(times));
motile = zeros(size(times));
nonmotile = zeros(size(times));
for curr_time = 1:length(times),
    for curr_cell = 1:length(all_data),
        temp = all_data{curr_cell};
        temp_times = [NaN; temp(2:end,22)];
        x = find(temp_times == times(curr_time));
        if ~isempty(x),
            if x == 2,
                cur_vel = temp(2, 19);
            else
                cur_vel = temp(x, 19);
            end
            if cur_vel > min_inst_vel,
                motile(curr_time) = motile(curr_time)+1;
            else
                nonmotile(curr_time) = nonmotile(curr_time)+1;
            end
            pop_vel{curr_time} = cat(1, pop_vel{curr_time}, cur_vel);
        end
        
    end
end
 
%change in cell motility over time
figure (cell_mot_chg_fig);
hold on;
plot(times, 100*(motile./(motile+nonmotile)),'b-')
title('Change in Population Cell Motility over Time');
xlabel('Time');
ylabel('Percent Motile Cells');
%create array to output
Motile_array = cat(2,times,(100*(motile./(motile+nonmotile))));
%output Motile array
%k = findstr(infile,'/');
%infile = infile(k(end)+1:end-12);
filename = strcat(infile,'motile_pop.txt');
dlmwrite(filename, Motile_array, '\t');
legend(strrep(filelist,'_','\_'));



%Construct and Plot Population inst. 2D velocity plots
pop_vel_mean = zeros(size(times));
pop_vel_error = zeros(size(times));
for curr_time = 1:length(times),
    pop_vel_mean(curr_time) = mean(pop_vel{curr_time});
    pop_vel_error(curr_time) = std(pop_vel{curr_time})/sqrt(length(pop_vel(curr_time)));
end
figure;
hold on;
plot(times, pop_vel_mean, 'b-');
plot(times, pop_vel_mean+pop_vel_error, 'r--');
plot(times, pop_vel_mean-pop_vel_error, 'r--');
title('Overall Population Mean +/- SEM Instantaneous 2D Velocity Over Time');
xlabel('Time');
ylabel('Instantaneous 2D Velocity');

%Plot the motile population's mean inst. 2D velocity
pop_vel_mean = zeros(size(times));
pop_vel_error = zeros(size(times));
for curr_time = 1:length(times),
    temp = pop_vel{curr_time};
    above_thresh_vel = temp(find(temp>=min_inst_vel));
    pop_vel_mean(curr_time) = mean(above_thresh_vel);
    pop_vel_error(curr_time) = std(above_thresh_vel)/sqrt(length(above_thresh_vel));
end
figure;
hold on;
plot(times, pop_vel_mean, 'b-');
plot(times, pop_vel_mean+pop_vel_error, 'r--');
plot(times, pop_vel_mean-pop_vel_error, 'r--');
title('Motile Population Mean +/- SEM Instantaneous 2D Velocity Over Time');
xlabel('Time');
ylabel('Instantaneous 2D Velocity');
%create array to output
TwoDIV_array = cat(2,times,pop_vel_mean,pop_vel_error);
%output 3DIV array
filename = strcat(infile,'2DIVmotpop.txt');
dlmwrite(filename, TwoDIV_array, '\t');


% Calculate and Plot 3D Inst Vel over Time

all_times = [];
for i = 1:length(all_data),
    temp = all_data{i};
    all_times = [all_times; temp(2:end,22)];    
end
times = unique(all_times(find(~isnan(all_times))));
pop_vel = cell(size(times));
motile = zeros(size(times));
nonmotile = zeros(size(times));
for curr_time = 1:length(times),
    for curr_cell = 1:length(all_data),
        temp = all_data{curr_cell};
        temp_times = [NaN; temp(2:end,22)];
        x = find(temp_times == times(curr_time));
        if ~isempty(x),
            if x == 2,
                cur_vel = temp(2, 10);
            else
                cur_vel = temp(x, 10);
            end
            if cur_vel > min_inst_vel,
                motile(curr_time) = motile(curr_time)+1;
            else
                nonmotile(curr_time) = nonmotile(curr_time)+1;
            end
            pop_vel{curr_time} = cat(1, pop_vel{curr_time}, cur_vel);
        end
        
    end
end


%Construct and Plot Population inst. 3D velocity plots
pop_vel_mean = zeros(size(times));
pop_vel_error = zeros(size(times));
for curr_time = 1:length(times),
    pop_vel_mean(curr_time) = mean(pop_vel{curr_time});
    pop_vel_error(curr_time) = std(pop_vel{curr_time})/sqrt(length(pop_vel(curr_time)));
end
figure;
hold on;
plot(times, pop_vel_mean, 'b-');
plot(times, pop_vel_mean+pop_vel_error, 'r--');
plot(times, pop_vel_mean-pop_vel_error, 'r--');
title('Overall Population Mean +/- SEM Instantaneous 3D Velocity Over Time');
xlabel('Time');
ylabel('Instantaneous 3D Velocity');
%create array to output
ThreeDIV_array = cat(2,times,pop_vel_mean,pop_vel_error);
%output ThreeDIV array
filename = strcat(infile,'3DIVpop.txt');
dlmwrite(filename, ThreeDIV_array, '\t');


%Plot the motile population's mean inst. 3D velocity
pop_vel_mean = zeros(size(times));
pop_vel_error = zeros(size(times));
for curr_time = 1:length(times),
    temp = pop_vel{curr_time};
    above_thresh_vel = temp(find(temp>=min_inst_vel));
    pop_vel_mean(curr_time) = mean(above_thresh_vel);
    pop_vel_error(curr_time) = std(above_thresh_vel)/sqrt(length(above_thresh_vel));
end
figure;
hold on;
plot(times, pop_vel_mean, 'b-');
plot(times, pop_vel_mean+pop_vel_error, 'r--');
plot(times, pop_vel_mean-pop_vel_error, 'r--');
title('Motile Population Mean +/- SEM Instantaneous 3D Velocity Over Time');
xlabel('Time');
ylabel('Instantaneous 3D Velocity');
%create array to output
ThreeDIVmot_array = cat(2,times,pop_vel_mean,pop_vel_error);
%output Motile ThreeDIV array
filename = strcat(infile,'3DIVmotpop.txt');
dlmwrite(filename, ThreeDIVmot_array, '\t');

warning on all;